% This code produces figures A3 and A4.

% parameters
a = 1;
b = 1; 
alpha=.45;
beta=.2;
H0 = 1;
lamda=0.5;

% chinese house demand
step=0.002;
C_chn = 0:step:0.06;

% intial condition
P_N_0 = b/a;
P_H_0 = (1-alpha-beta)/(alpha+beta)*b/H0;
w = b;
e_N_0 = alpha/(alpha+beta);
e_T_0 = beta/(alpha+beta);

% commute cost distribution
phi_mu=0.032;phi_sd=0.008; %normal distribution


% equilibrium
soln=zeros(1,size(C_chn,2));
for i=1:size(C_chn,2)
chn=C_chn(i);
syms x
e=normcdf(x,phi_mu,phi_sd);
f=1-((1-chn/((alpha+beta)*H0))*...
    ((1+e)/(1-e))...
    )^(1-alpha-beta)-x ;
soln(i)=vpasolve(f,x); 
end

phi_bar=soln;

e_bar = normcdf(phi_bar,phi_mu,phi_sd); 

P1_H = (1-alpha-beta)*b*(1-e_bar)./((alpha+beta)*H0-C_chn);
P0_H = (1-alpha-beta)*b*(1+e_bar)./((alpha+beta)*H0);
e1_N = alpha/b*P1_H*H0+alpha*(1-e_bar*lamda);
e0_N = alpha/b*P0_H*H0+alpha*(1+e_bar*lamda);
e1_T = 1+e_bar-e1_N;
e0_T = 1-e_bar-e0_N;
e1 = 1+e_bar;
e0 = 1-e_bar;

% normalize variables
P1_H = P1_H./P1_H(1);
P0_H = P0_H./P0_H(1);
e1_N = e1_N./e1_N(1);
e0_N = e0_N./e0_N(1);
e1_T = e1_T./e1_T(1);
e0_T = e0_T./e0_T(1);
e1 = e1./e1(1);
e0 = e0_T./e0(1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure a4
subplot(2,3,1)
plot(C_chn,e1,C_chn,e0, 'r--','LineWidth',1)
legend('Region 1','Region 0','Location','southwest')
xlabel('Foreign demand ($C_f^H$)','Interpreter','latex')
ylabel('$e$ (normalized by initial equilibrium)','Interpreter','latex')
title('Total Employment ($e$)','Interpreter','latex')

subplot(2,3,2)
plot(C_chn,e1_N,C_chn,e0_N, 'r--','LineWidth',1)
ylim([1 1.02])
xlabel('Foreign demand ($C_f^H$)','Interpreter','latex')
ylabel('$e^N$(normalized by initial equilibrium)','Interpreter','latex')
title('Nontradable-sector Employment ($e^N$)','Interpreter','latex')

subplot(2,3,3)
plot(C_chn,e1_T,C_chn,e0_T, 'r--','LineWidth',1)
xlabel('Foreign demand ($C_f^H$)','Interpreter','latex')
ylabel('$e^T$(normalized by initial equilibrium)','Interpreter','latex')
title('Tradable-sector Employment ($e^T$)','Interpreter','latex')

subplot(2,3,4)
plot(C_chn,P1_H,C_chn,P0_H, 'r--','LineWidth',1)
xlabel('Foreign demand ($C_f^H$)','Interpreter','latex')
ylabel('$P^H$(normalized by initial equilibrium)','Interpreter','latex')
title('Housing Price ($P^H$)','Interpreter','latex')

subplot(2,3,5)
plot(C_chn,e_bar,'LineWidth',1)
ylim([0 0.03])
xlabel('Foreign demand ($C_f^H$)','Interpreter','latex')
ylabel('$\delta$','Interpreter','latex')
title('Commuters ($\delta$)','Interpreter','latex')

print(gcf,'..\results\figureA4.png','-dpng','-r300'); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% different phi_mu

% baseline mu
e1_mu03=e1;
e1_N_mu03=e1_N;
e1_T_mu03=e1_T;
P1_H_mu03=P1_H;
e_bar_mu03=e_bar;

% mu=0.05
phi_mu=0.05;phi_sd=0.008; %normal distribution

% equilibrium
soln=zeros(1,size(C_chn,2));
for i=1:size(C_chn,2)
chn=C_chn(i);
syms x
e=normcdf(x,phi_mu,phi_sd);
f=1-((1-chn/((alpha+beta)*H0))*...
    ((1+e)/(1-e))...
    )^(1-alpha-beta)-x ;
soln(i)=vpasolve(f,x); 
end

phi_bar=soln;

e_bar = normcdf(phi_bar,phi_mu,phi_sd); 

P1_H = (1-alpha-beta)*b*(1-e_bar)./((alpha+beta)*H0-C_chn);
P0_H = (1-alpha-beta)*b*(1+e_bar)./((alpha+beta)*H0);
e1_N = alpha/b*P1_H*H0+alpha*(1-e_bar*lamda);
e0_N = alpha/b*P0_H*H0+alpha*(1+e_bar*lamda);
e1_T = 1+e_bar-e1_N;
e0_T = 1-e_bar-e0_N;
e1 = 1+e_bar;
e0 = 1-e_bar;

% normalize variables
P1_H = P1_H./P1_H(1);
P0_H = P0_H./P0_H(1);
e1_N = e1_N./e1_N(1);
e0_N = e0_N./e0_N(1);
e1_T = e1_T./e1_T(1);
e0_T = e0_T./e0_T(1);
e1 = e1./e1(1);
e0 = e0_T./e0(1);

e1_mu05=e1;
e1_N_mu05=e1_N;
e1_T_mu05=e1_T;
P1_H_mu05=P1_H;
e_bar_mu05=e_bar;


%mu=0.02
% commute cost distribution
phi_mu=0.02;phi_sd=0.008; %normal distribution

% equilibrium
soln=zeros(1,size(C_chn,2));
for i=1:size(C_chn,2)
chn=C_chn(i);
syms x
e=normcdf(x,phi_mu,phi_sd);
f=1-((1-chn/((alpha+beta)*H0))*...
    ((1+e)/(1-e))...
    )^(1-alpha-beta)-x ;
soln(i)=vpasolve(f,x); 
end

phi_bar=soln;

e_bar = normcdf(phi_bar,phi_mu,phi_sd); 

P1_H = (1-alpha-beta)*b*(1-e_bar)./((alpha+beta)*H0-C_chn);
P0_H = (1-alpha-beta)*b*(1+e_bar)./((alpha+beta)*H0);
e1_N = alpha/b*P1_H*H0+alpha*(1-e_bar*lamda);
e0_N = alpha/b*P0_H*H0+alpha*(1+e_bar*lamda);
e1_T = 1+e_bar-e1_N;
e0_T = 1-e_bar-e0_N;
e1 = 1+e_bar;
e0 = 1-e_bar;

% normalize variables
P1_H = P1_H./P1_H(1);
P0_H = P0_H./P0_H(1);
e1_N = e1_N./e1_N(1);
e0_N = e0_N./e0_N(1);
e1_T = e1_T./e1_T(1);
e0_T = e0_T./e0_T(1);
e1 = e1./e1(1);
e0 = e0_T./e0(1);

e1_mu02=e1;
e1_N_mu02=e1_N;
e1_T_mu02=e1_T;
P1_H_mu02=P1_H;
e_bar_mu02=e_bar;

% figure a5
subplot(2,3,1)
plot(C_chn,e1_mu05,'black-.',C_chn,e1_mu03,'blue-',C_chn,e1_mu02,'m--','LineWidth',1)
legend('\mu=0.05','\mu=0.032','\mu=0.02','Location','northwest')
xlabel('Foreign demand ($C_f^H$)','Interpreter','latex')
ylabel('$e$ (normalized by initial equilibrium)','Interpreter','latex')
title('Total Employment ($e$)','Interpreter','latex')

subplot(2,3,2)
plot(C_chn,e1_N_mu05,'black-.',C_chn,e1_N_mu03,'blue-',C_chn,e1_N_mu02, 'm--','LineWidth',1)
xlabel('Foreign demand ($C_f^H$)','Interpreter','latex')
ylabel('$e^N$(normalized by initial equilibrium)','Interpreter','latex')
title('Nontradable-sector Employment ($e^N$)','Interpreter','latex')

subplot(2,3,3)
plot(C_chn,e1_T_mu05,'black-.',C_chn,e1_T_mu03,'blue-',C_chn,e1_T_mu02, 'm--','LineWidth',1)
xlabel('Foreign demand ($C_f^H$)','Interpreter','latex')
ylabel('$e^T$(normalized by initial equilibrium)','Interpreter','latex')
title('Tradable-sector Employment ($e^T$)','Interpreter','latex')

subplot(2,3,4)
plot(C_chn,P1_H_mu05,'black-.',C_chn,P1_H_mu03,'blue-',C_chn,P1_H_mu02, 'm--','LineWidth',1)
xlabel('Foreign demand ($C_f^H$)','Interpreter','latex')
ylabel('$P^H$(normalized by initial equilibrium)','Interpreter','latex')
title('Housing Price ($P^H$)','Interpreter','latex')

subplot(2,3,5)
plot(C_chn,e_bar_mu05,'black-.',C_chn,e_bar_mu03,'blue-',C_chn,e_bar_mu02,'m--','LineWidth',1)
xlabel('Foreign demand ($C_f^H$)','Interpreter','latex')
ylabel('$\delta$','Interpreter','latex')
title('Commuters ($\delta$)','Interpreter','latex')

print(gcf,'..\results\figureA5.png','-dpng','-r300'); 

